Phases of attractive spin-imbalanced fermions in square optical lattices 
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We determine the relative stability of different ground-state phases of spin-imbalanced popula- 
tions of attractive fermions in square lattices. The phases are systematically characterized by the 
symmetry of the order parameter and the real- and momentum-space structures using Hartree- 
Fock-Bogoliubov theory. We find several type of unidirectional Larkin-Ovchinikov-type phases. We 
discuss the effect of commensuration between the ordering wave vector and the density imbalance, 
and describe the mechanism of Fermi surface reconstruction and pairing for various orders. A robust 
supersolid phase is shown to exist when the ordering wave vector is diagonally directed. 
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There has been a surge of interest in the possibility of 
realizing unconventional fermionic superfluids using cold 
atomic gases. Amongst the many possibilities offered 
by the highly tunable Hamiltonians available in cold- 
atoms experiments, the simplest remains that of a two- 
component density-imbalanced populations of fermions 
with attractive interactions. The theoretical study of 
such systems dates back to Fulde and Ferrel (FF) [Tj and 
Larkin and Ovchinikov (LO) [2] who independently sug- 
gested that the mismatch between the Fermi surfaces of 
the two species could result in the formation of a conden- 
sate of finite- momentum pairs. Atomic gases offer a di- 
rect route to the realization of FFLO phases, circumvent- 
ing most of the difficulties of solid state systems, thanks 
to the possibility of controlling independently the den- 
sity of the two species, the absence of disorder and, most 
importantly, the ability to engineer strong interactions. 
In spite of this, the existence of an FFLO phase in two 
and three dimensions has been argued to be confined to a 
small range of interaction strengths and polarizations 0, 
and detection has remained elusive. 

The possibility of using optical lattices has been sug- 
gested by several authors |S1 15] as a key ingredient to ob- 
serve FFLO-type states. The best empirical indication 
that this may be the case is provided by experiments on 
strongly-correlated-electrons materials and the fact that, 
when doped, these systems show a tendency toward for- 
mation of inhomogeneities in the form of spin, charge 
and, possibly, pairing density waves. The relevance of 
these experiments to the properties of attractive fermions 
in optical lattices stems from the belief that, in both 
cases, the essential physics can be captured by a one- 
band Hubbard model: with an on-site repulsive inter- 
action for many of the electronic systems and an on-site 
attraction in an optical lattice. The attractive and repul- 
sive cases are mapped into each other by a particle-hole 
transformation [4 and the presence of spin-texture in the 
doped repulsive case translates into the occurrence of a 
modulated superfiuid in an imbalanced population of at- 
tractive fermions, ie. an FFLO phase. This is reenforced 
by recent quantum Monte Carlo results [6] on the two- 
dimensional repulsive model showing spin-density waves 



with long wavelength modulation. 

Despite this mapping and several works addressing 
the existence of a possible FFLO phase [7-^, the na- 
ture of the ground state phases in a spin-imbalanced 
two-dimensional optical lattice remains largely undeter- 
mined. On the one hand, information on the repul- 
sive model is entirely confined to the case of unpolar- 
ized systems, which maps into the attractive case at 
half- filling: + = 1; for the case of imbalanced 
fermionic population, one is interested in the more gen- 
eral case of a polarized system and arbitrary density. On 
the other hand, works addressing the physics in the lat- 
tice have either focused on the single plane-wave form of 
the order parameter [71 [8] , which is likely incorrect at low 
temperature [2], or on selected states [9 {e.g. fixed modu- 
lation wavelength) because of the challenge of removing 
large finite-size effects. 

In this paper, we determine systematically the ground 
state phases of the attractive two-dimensional spin- 
imbalanced optical lattice system away from half-filling. 
Within mean-field theory, we converge the numerical cal- 
culations to the thermodynamic limit, and character- 
ize the order parameter and the real- and momentum- 
space properties for different densities and polarizations. 
Small to moderate interaction strengths are considered, 
where mean-field theory is expected to capture the cor- 
rect physics. We find the existence of several types of 
LO phases in a large region of the parameter space, with 
distinct physical characteristics which have not been de- 
scribed before. Particularly noteworthy is a robust su- 
persolid phase obtained when the ordering wave vector 
is directed along the diagonal direction. 

Results presented in this work are obtained using 
Hartree-Fock-Bogoliubov theory so that modulations in 
charge, spin and pairing are all handled on the same foot- 
ing. The starting Hamiltonian reads 

H = -t^ c\^Cj^ Ui^rii^ - jarii - ^rrii, (1) 

{ij)a i 

where are fermionic annihilation operators of spin cr 
on site i, n^^ = 4^Q^, = n^^ + n^;, = n^^ - n^;. 
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In order to accommodate the inhomogeneities, the cal- 
culations are performed on finite simulation cells whose 
shape is dictated by the symmetry of the targeted phase. 
The cells are characterized by two basis vectors, Li 
and L2, whose components are integers. Once the cell 
shape is chosen, we define Bloch states as Cj{k) cx 
Cj+L exp \ik ' L\ where L = niLi + n2l/2, and A: is a 
vector that varies freely within the first BZ of the recip- 
rocal lattice. Then, using these states and the mean-field 
approximation, the Hamiltonian decouples into a sum of 
/c-dependent pieces, H = H{k)^ of the form 



H{k) = [c|cj 
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where and represent an array (row) of operators, 
{ci^(/c)} and — k)} respectively, with index i run- 

ning over the TV sites of the simulation cell, and G is 
defined below. H and A are N x N matrices with ele- 
ments 



-tij{k) + 5ij{Di, 
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In the above, tij{k) = '^j^ex.p{ik • L)tij-^L^ s 
and Dicj^ A^, fi and h are determined by the requirement 
that the Free energy F = (H) — TS is a minimum for 
the target average densities Ucr. This amounts to the 
following self-consistency equations 



A-. =U J dk{cl{k)ci,{k)) , 
A, = -U [ dk{cl{k)cl_M), 



(4) 



-N-' ^ I dk{cl{k)c,,{k)) ^ 



with expectation values evaluated in the simulation cell. 

The variable G in Eq. Q is a vector such that = G-L 
gives the twist angle of the pairing order parameter un- 
der translation by L. For collinear phases, is either 
or TT, giving periodic or anti-periodic boundary condi- 
tions on (ct^cj^). Note that uniform phases with a spiral 
order parameter (FF type) amount to the choice of a one- 
site cell and a G equal to the wave- vector of the spiral 
modulation. Thus H{k) is specified by a 2 x 2 matrix, 
and the eigenvalues and eigenvectors can be computed 
analytically, in the spiral phase [1. 

Given a set of symmetry-equivalent q vectors, linear 
response theory can be used to show that the onset of 
instabilities of the form 



(5) 



must happen at exactly the same value of [/, regardless 
of the choice of A^. In order to determine the correct 




FIG. 1: Left Panel: local order parameter, max^ | (c^^Ci^) |, 
and leading wavevector (inset) versus U for (ci^Cii) oc e**^"^'^^ 
(Spiral), (ci^Cii) oc cos{qx • n) (Linear), (ci^Ci^) oc cos{qx • 
ri)-\-cos{qyri) (Checkbd) with = |^|(1,0) and qy = |^|(0,1). 
Right panel: relative energies (units of t) of the three phases. 



form of order parameter, we proceed as follows. We first 
determine Uc and the associated non-zero wave-vector 
qc using the single plane-wave form as this allows for 
a quick exploration of phase space. We find that qc is 
directed along any of the four, symmetry-equivalent di- 
rections (±1,0), (0, ±1) and, therefore, any linear com- 
bination of the four associated plane waves is a candi- 
date ground state order parameter just above Uc. To 
resolve which one leads to the largest lowering of energy, 
we proceed by solving the mean-field equations for the 
three cases corresponding to spiral (A^ oc exp{iqx ■ r^)), 
unidirectional (A^ cx cos{qx-ri)) and checkerboard (A^ cx 
cos{qx-ri)-\-cos{qy-ri)) pairing density wave and track the 
evolution of |g| as a function of U. Explicit calculations in 
large simulations cells on the repulsive model (with sim- 
ulated annealing starting with random initial fields [10^) 
have shown that instabilities involving vectors in differ- 
ent, non-equivalent directions, which the above approach 
would miss, are unlikely to occur in the range of U con- 
sidered here. 

Mean- field results [10 on the repulsive model and the 
particle-hole transformation relating the attractive and 
repulsive models imply the existence of the following 
properties at half- filling, i.e., when the average parti- 
cle density is precisely one fermion per site: 1) a crit- 
ical U exists such that, above it, the system develops 
a phase with an inhomogeneous order parameter 2) the 
pairing order parameter is characterized by a wave vector 
\q\ = mir where m = — is the average magnetiza- 
tion 3) both fermionic species have a gap in their single 
particle spectrum 4) as U grows larger there is a transi- 
tion from a pair-density- wave in the (1,0) state to one in 
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FIG. 2: Local properties and gradient of the momentum distribution for — n^ — —0.05 at two densities: (A) n — 0.96 (top 
row) and (B) n — 0.93 (bottom). In both cases, the finite-momentum of the pair, is in the (1, 0)-direction. The inset in the 
first column shows the evolution of A = rmi/q as a function of density; case A belongs to the commensurate regime with unit 
density of excess-spin particles per node, while case B is in the incommensurate regime. The third column is the local density 
of state measured at site 0. In the fourth column, each panel shows two halves of the Fermi surface, for the t (minority) and 
\, (majority) spins, respectively. The arrows in the figure indicate the pairing construction and applies to every pair across the 
FS, leading to one common q. 



the (1,1) direction 5) order can be arbitrarily broken into 
a charge or a pairing instability or a combination of the 
two. This last point is a consequence of the possibility 
of breaking spin symmetry in any of the three equivalent 
direction in the repulsive case. It implies that charge and 
pairing orders compete, in the sense that the larger one 
is, the smaller the other must be. 

We next present results away from half- filling. Figure [l] 
shows that the ground state has unidirectional LO order, 
with illustrative examples for m = — = —0.095 
and n = 0.95. The presence of the lattice is manifested 
as the significantly different shape of the non-interacting 
Fermi surfaces when compared to their circular counter- 
parts in the continuum. Such differences are largest when 
Tier = 0.5 and gradually disappear as the density is re- 
duced. For this reason, we repeat a similar analysis in dif- 
ferent regimes of density and polarization, n = 0.95, 0.75 
and 0.5 and m/n = 0.1 and 0.4, so as to have a rather 
complete picture close and away from half-filling, at small 
and large polarization and in an interaction range that 
extends from Uc up to U = 4. The results are summa- 
rized in Table [H These are consistent with earlier results 
addressing the physics of FFLO phases in the continuum, 
which found checkerboard order close to Uc thanks to an 
expansion in powers of the order parameter [11, but show 
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TABLE I: Critical parameter at the onset of order for a few 
densities (n) and polarizations (m/n). The wavevector is di- 
rected along (1,0). The two letters in last line in each table 
cell describe the type of order just above Uc and at [/ = 4, 
respectively (L=linear, C=checkerboard). 



that for n > 0.75 lattice effects are strong enough to re- 
turn unidirectional order independently of polarization 
or proximity to Uc Given the difficulty in observing un- 
conventional pairing in the continuum, the latter is the 
density regime where an experimental realization of the 
FFLO state should be most feasible. We will therefore 
focus on it in the following. 

Figure [2] characterizes the FFLO phase at = 3t 
and for m = 0.05 in two qualitatively different density 
regimes. In particular, for n > 0.95, the wave vector is 
precisely determined by the magnetization via the same 
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relation holding at half- filling, q = mir. This commen- 
surate regime is characterized by a density of one excess 
particle per node (of the order parameter) and conse- 
quently band-insulating behavior along the node. This is 
most clearly seen in the local density of states at the node, 
with both species showing a gap at the Fermi energy. 
Correspondingly the gradient of the momentum distri- 
bution shows no sharp lines indicative of the existence 
of a Fermi surface. The commensurate regime here has, 
however, some important distinctions from half-filling. 
First, the interchangeability between pairing and charge 
orders is broken as soon as the average density deviates 
from n = 1, and pairing emerges as the dominant order. 
Second, the density is not perfectly uniform, as it is for 
the purely superfiuid phase at half-filling. The density is 
instead characterized by a weak modulation that reflects 
the different degrees of localization of the Andreev states 
for the two spin species. The density profile shows weak 
peaks at the nodes indicating that the majority species 
has stronger localization. 

In the second density regime, n < 0.95, the major- 
ity spin species develop a finite density of states at the 
nodes of the order parameter. Fermi arcs (sharp lines) 
are seen in |Vn/c||. The arc in Fig.[2]B is of perfectly one- 
dimensional nature, indicating a complete decoupling be- 
tween the metallic states living at different nodes. At 
larger polarization, the arcs will more strongly resem- 
ble the underlying Fermi surface of the non-interacting 
species. Finally, the density profile shows minima at 
the nodes, rather than the maxima observed in case A, 
as a direct consequence of the gradual emptying of the 
majority-spin Andreev bands visible in the local density 
of states. This is a potentially important experimental 
signature as it allows one to characterize the nature of 
the nodal metal via a static local property instead of a 
collective property such as the distance between nodes or 
the presence of Fermi arcs. 

Next we consider the phases at larger interaction 
strengths U. Results on the repulsive model [10 indi- 
cate that the ordering wave vector switches to the (1,1) 
direction at sufficiently small m and U > 3t. We focus on 
U = At for an extensive and systematic examination of 
the properties of these phases. Indeed unidirectional LO 
states along the diagonal direction are found; however, we 
find that the system can accommodate both pairing and 
charge orders as non- competing instabilities away from 
half-filling. The real- and momentum-space properties 
of the phases are shown in Fig. |3] In contrast with the 
(l,0)-direction paired state, pairing here is achieved by a 
reconstruction of the Fermi surfae. The non-interacting 
surface of the majority spin is stretched, while that of 
the minority spin is compressed, along the (-1,1) direc- 
tion so that the reconstructed surfaces become "nested" 
via /c ^ — /c + g as illustrated in the middle row of Fig|3] 
The charge-density wave that develops at larger density 
is a consequence of the (7r,7r) nesting along the (1,1)- 




FIG. 3: Left panels: Evolution of charge and pairing order as 
density is reduced, for m = 0.05. Right panels: corresponding 
evolution of the gradient of the momentum distributions, and 
illustration of Fermi surface reconstruction and pairing. The 
pairing wave-vector, g, is in the ( — 1, l)-direction. The solid 
(green) lines give the non-interacting Fermi surface. 

direction which results from the reconstruction; its am- 
plitude increases with the amount of nested states. The 
charge order is complementary to pairing, and further 
lowers the energy. 

The new phase discussed above has an important dis- 
tinction from the supersolid phase that can occur at half- 
filling. There both charge and pairing orders are associ- 
ated with wave- vectors coupling the same regions around 
the Fermi surface: if pairing order happens at charge 
order must happen at (tt, tt) — g and one can dial the form 
of the order parameter interpolating between a purely 
superfiuid and a purely charge-density wave state. Here, 
away from half- filling, however, pairing is characterized 
by a polarization dependent wave-vector q while charge 
order appears at (tt, tt) and cannot be dialed away. 

Our results are summarized in the phase diagram of 
Fig. |4] for /7 = 3t and U = 4t. We found no evidence 
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FIG. 4: Uji = 3 (left) and U jt = 4 (right) T = phase dia- 
gram. The color axis reports values of A = mi^ l\q\. Triangles 
represent the commensurate phases (A = 1). Up-triangles 
have q oc (1,1) (DC=Diagonal commensurate) while down- 
triangles have ^ oc (1,0) (VC= Vertical Commensurate). Cir- 
cles and square are Vertical and Diagonal Incommensurate 
phases (VI and DI). 



of diagonal order at = 3t. Its region of stability for 
the commensurate phase is limited to small polarization 
and density close to one. At = 4t, the diagonal wave 
vector is almost always commensurate with m, and pair- 
ing and charge order coexist in all but a small fraction of 
the phase diagram where the diagonal phase is the cor- 
rect ground state. This suggests that the charge-density 
wave play a small but important role in stabilizing di- 
agonal order. Although the different symmetry between 
diagonal and vertical phases implies that the transition 
is discontinuous and accompanied by phase separation, 
we have not attempted an in-depth study of how this af- 
fects the situation at = 4t. However, a simpler analysis 
based on considering phase separation in two phases hav- 
ing the same density and different magnetization leads to 
a narrow coexistence region (Am = 0.02 at n = 1) and 
makes it sensible to assume that the broad feature of the 
phase diagram are, in fact, robust. 

In summary, we have determined the type of 
phases that an imbalanced population of two fermionic 
species with attractive interactions support in the two- 
dimensional square lattice. Their real- and momentum- 
space properties are quantitatively characterized. We 
find that unidirectional LO states are the most sta- 
ble mean field solutions in a large range of parameter 
regimes, and checkerboard order is only clearly favored 
at small density and large polarization. At lower [/, the 
finite- momentum pairing wave- vector is along (1,0). We 



have shown the insulating versus conducting nature of 
the nodal region of the superconducting order parameter, 
and its interplay with commensuration effects. Related 
to these effects, a new phase with supersolid order is seen 
when the ordering wave vector is directed along the (1,1) 
direction. Our results suggest that, besides time of flight 
or Bragg spectroscopy, the different phases we discussed 
are also identifiable by their local density profiles. This, 
in turn, should help their real-space characterization in 
the presence of a confining potential. 
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